﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

//using MathWorks.MATLAB.NET.Arrays;
//using InverseMatrix;
//using Det;
//using DivideMatrix;

using DotNumerics.LinearAlgebra;

namespace HAMath
{  
    public static class MathHleper
    {

        public static double GetDistance(KrigingPoint pt1, KrigingPoint pt2)
        {
            return Math.Sqrt((pt2.X - pt1.X) * (pt2.X - pt1.X) + (pt2.Y - pt1.Y) * (pt2.Y - pt1.Y));
        }


        /// <summary>
        /// 判断一个排列是奇排列还是偶排列。是偶排列返回true。没有用着
        /// </summary>
        /// <param name="values"></param>
        /// <returns></returns>
        private static bool IsEvenOrder(List<int> values)
        {
            int k = 0;

            for (int i = 0; i < values.Count - 1; i++)
            {
                for (int j = i + 1; j < values.Count; j++)
                {
                    if (values[j] < values[i])
                    {
                        k++;
                    }
                }
            }

            int residue; //余数
            Math.DivRem(k, 2, out residue);

            if (residue == 0)
            {
                return true;
            }

            return false;

        }


        /// <summary>
        /// 求一个方阵的行列式的元素。没有用着
        /// </summary>
        /// <param name="data"></param>
        /// <returns></returns>
        private static double[] DeterminantElement(double[,] data)
        {
            //if (data.GetLength(0) != data.GetLength(1)) return double.NaN;
            //if (data.GetLength(0) < 2) return double.NaN;

            int num = data.GetLength(0);
            double[] result;

            if (num == 2)
            {
                result = new double[2];
                result[0] = data[0, 0] * data[1, 1];
                result[1] = data[0, 1] * data[1, 0];
                return result;
            }

            int length = 1;
            int temNum = num;
            while (temNum > 1)
            {
                length = length * temNum;
                temNum--;
            }
            result = new double[length];

            for (int i = 0; i < num; i++)
            {
                double[,] temData = new double[num - 1, num - 1];

                for (int m = 1; m < num; m++)
                {
                    for (int n = i + 1; n < i + num; n++)
                    {
                        int index1 = m - 1, index2 = n - i - 1;
                        int index3 = m, index4 = n;
                        if (index4 >= num) index4 = index4 - num;

                        temData[index1, index2] = data[index3, index4];
                    }
                }

                double[] temResult = DeterminantElement(temData);

                for (int j = 0; j < temResult.Length; j++)
                {
                    result[temResult.Length * i + j] = data[0, i] * temResult[j];
                }

            }

            return result;

        }


        /// <summary>
        /// 求矩阵的行列式。矩阵必须是方形矩阵
        /// </summary>
        /// <param name="data"></param>
        /// <returns></returns>
        private static double GetDeterminant(double[,] a, int k)
        {
            double result = 0;

            int c, d;
            double temp;
            int i, j, x, y;
            double[,] temp1;
            double[,] temp2;

            temp1 = new double[k, k];

            for (int m = 0; m < k; m++)
            {
                for (int n = 0; n < k; n++)
                {
                    temp1[m, n] = a[m, n];
                }
            }

            d = k;

            for (c = 0; c < d; c++)
            {
                a = new double[d, d]; // '在每降一阶的循环中，a()初始化为原传入的数组
                for (i = 0; i < d; i++) //\\ To d
                {
                    for (j = 0; j < d; j++) // To d
                    {
                        a[i, j] = temp1[i, j];
                    }
                }

                k = d; //'初始化k,每调用一次ConvertA,k值减1

                if (k >= 2)
                {

                    //''求行列式a()的第一行元素a(1,c)的余子式

                    temp2 = new double[k - 1, k - 1]; //'求行列式a()的第一行元素a(1,c)的余子式

                    for (x = 0; x < k - 1; x++)  //To k - 1
                    {
                        for (y = 0; y < k - 1; y++) // To k - 1
                        {
                            if (y < c) // Then
                            {
                                temp2[x, y] = a[x + 1, y];
                            }
                            else
                            {
                                temp2[x, y] = a[x + 1, y + 1];
                            }
                        }
                    }

                    k = k - 1;

                    a = new double[k, k];
                    for (x = 0; x < k; x++) //To k
                    {
                        for (y = 0; y < k; y++)  // To k
                        {
                            a[x, y] = temp2[x, y];
                        }
                    }
                    //''以上为求行列式a()的第一行元素a(1,c)的余子式


                    temp = Math.Pow(-1, c) * temp1[0, c] * GetDeterminant(a, k); ; // 'k为阶数，c为元素第几列

                    result = result + temp;
                }
                else if (k == 1)  // Then
                {
                    result = result + a[0, 0];
                }


            }

            return result;

        }


        /// <summary>
        /// 求参数矩阵的伴随矩阵
        /// </summary>
        /// <param name="data"></param>
        /// <returns></returns>
        private static double[,] GetAdjointMatrix(double[,] matrix)
        {
            if (matrix.GetLength(0) != matrix.GetLength(1)) return null;
            if (matrix.GetLength(0) < 2) return null;

            double[,] result = new double[matrix.GetLength(0), matrix.GetLength(1)];

            int num = matrix.GetLength(0);

            int k = 0;  //计数器

            for (int i = 0; i < num; i++)
            {
                for (int j = 0; j < num; j++)
                {
                    double[,] temMatrix = new double[num - 1, num - 1];
                    for (int m = 0; m < num; m++)
                    {
                        if (m == i) continue;
                        for (int n = 0; n < num; n++)
                        {
                            if (n == j) continue;

                            int index1 = m, index2 = n;

                            if (m > i) index1 = m - 1;
                            if (n > j) index2 = n - 1;

                            temMatrix[index1, index2] = matrix[m, n];

                        }
                    }

                    int residue; //余数
                    Math.DivRem(k, 2, out residue);

                    //if (residue == 0)
                    //{
                    double p = i + 1 + j + 1;
                    double t = Math.Pow(-1.0, p);
                    result[j, i] = GetDeterminant(temMatrix, temMatrix.GetLength(0)) * t;
                    //}
                    //else
                    //{
                    //    result[i, j] = -this.GetDeterminant(temData);
                    //}

                    k++;

                }
            }

            return result;

        }



        //private static InvMatrix invMatrix;
        //private static DetClass det;
        //private static DivideMartrix divMatrix;

        /// <summary>
        /// 求参数矩阵的逆矩阵
        /// </summary>
        /// <param name="data"></param>
        /// <returns></returns>
        public static double[,] GetInverseMatrix(double[,] matrix)
        {
            if (matrix.GetLength(0) != matrix.GetLength(1)) return null;
            if (matrix.GetLength(0) < 2) return null;


            DotNumerics.LinearAlgebra.Matrix dotMatrix = new Matrix(matrix);
            DotNumerics.LinearAlgebra.Matrix dotResult = dotMatrix.Inverse();

            
            double[,] result = dotResult.CopyToArray();

            return result;

            //下面的是matlab的方法

            //double[,] result = new double[matrix.GetLength(0), matrix.GetLength(1)];
                       

            //if (invMatrix == null)
            //{
            //    invMatrix = new InvMatrix();
            //}
            //if (det == null)
            //{
            //    det = new DetClass();
            //}
            
            
            //MWNumericArray m = matrix;            
            

            //MWNumericArray mResult = invMatrix.InverseMatrix(m) as MWNumericArray;

            //if (mResult.IsNumericArray)
            //{
            //  Array a=  mResult.ToArray(MWArrayComponent.Real);
            //  result = a as double[,];
            //}

            //m.Dispose();
            //mResult.Dispose();
            
            //return result;

        }
        

        /// <summary>
        /// 矩阵相乘
        /// </summary>
        /// <param name="data1"></param>
        /// <param name="data2"></param>
        /// <returns></returns>
        public static double[] MultiplyMatrix(double[,] matrix1, double[] matrix2)
        {
            if (matrix1.GetLength(0) != matrix1.GetLength(1)) return null;
            if (matrix1.GetLength(1) != matrix2.Length) return null;

            double[] result = new double[matrix2.Length];

            for (int i = 0; i < matrix1.GetLength(0); i++)
            {
                for (int j = 0; j < matrix1.GetLength(1); j++)
                {
                    result[i] = result[i] + matrix1[i, j] * matrix2[j];
                }
            }

            return result;
        }


        public static double[] DivideMatrix(double[] matrix1, double[,] matrix2)
        {
            //这是Matlab的方法，不用了

            //if (divMatrix == null)
            //{
            //    divMatrix = new DivideMartrix();
            //}

            //double[] result = new double[matrix1.Length];

            //MWNumericArray m1=matrix1 ;
            //MWNumericArray m2=matrix2 ;

            //MWNumericArray mResult = divMatrix.DivMatrix(m1, m2) as MWNumericArray;

            //if (mResult.IsNumericArray)
            //{
            //    Array a = mResult.ToArray(MWArrayComponent.Real);
            //    double[,] temResult = a as double[,];

            //    for (int i = 0; i < temResult.GetLength(1); i++)
            //    {
            //        result[i] = temResult[0, i];
            //    }
            //}

            //return result;

            return null;
        }


    }
}
